rm(list=ls())

load.module("glm")
load.module("lecuyer")


# Purpose: 
#  Run the JAGS model: MAIN MODEL. 
#
#
#
#



# LOAD  


load("./models/m301/for_jags.Rdata")
cmpemp <- for_jags$cmpemp
phi.select <- for_jags$phi.select 
theta.select <- for_jags$theta.select
rho.select <- for_jags$rho.select

###EC: this is the ebprior 

load("cee_prior.Rdata")
for_jags$mu0 <- ebprior$mu0
for_jags$sigmainv0 <- ebprior$sigmainv0

# DATA / CONSTRAINTS

mrg.y <- as.matrix(cmpemp[,c("lr1", "lr2", "lr3", "lr4", "lr5", "lr6", "lr7", "lr8")] )
mrg.y <- apply(mrg.y, 2, scale, scale=F)

P0 <- R0 <- T0 <- 0.01	
p0 <- r0 <- t0 <- 0
l0 <- 0
L0 <- 0.01
e0 <- f0 <- 0.001
b0 <- 0
B0 <- 1

sigmainv0 <- as.vector(for_jags$sigmainv0)

theta <- rep(NA, nrow(theta.select))
theta[theta.select$start==1] <- 0

rho <- rep(NA, nrow(rho.select))
rho[rho.select$entity=="EP"] <- 0

phi <- matrix(NA,nrow(phi.select),1)
phi <- as.vector(rbind(phi,for_jags$mu0))



# INITS

inits <- list( 
		list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=1234) ,
		list(".RNG.name"="lecuyer::RngStream", ".RNG.seed"=99999) 
	)


# ESTIMATE 
##########

jags <- jags.model('./models/mcss.bug',  inits=inits, data = list(
														'mrg.y' = mrg.y,
														'N' = nrow(mrg.y), 'L' = ncol(mrg.y),
														'mrg.it' = cmpemp$it,
														'mrg.ce' = cmpemp$ce,
														'mrg.c' = cmpemp$c,
														'phi.z' = phi.select$z, 
														'phi.w' = phi.select$w,
														'phi.IT' = nrow(phi.select), 
														'phi.i' = phi.select$i, 
														'phi.t' = phi.select$t,
														'phi.T' = phi.select$T + 1,
														'phi' = phi,
														'sigmainv0' = sigmainv0,
														'l0' = l0, 'L0' = L0, 
														'e0' = e0, 'f0' = f0,
														't0' = t0, 'T0' = T0,
														'r0' = r0, 'R0' = R0,
														'b0' = b0, 'B0' = B0,
														'mrg.CE' = length(unique(cmpemp$ce)), 
														'mrg.C' = length(unique(cmpemp$c)), 
														'mrg.I' = length(unique(phi.select$i)), 
														'rho' = rho,
														'theta' = theta
															 ), n.chains = 2)

update(jags, 10000)
phimcmc <- coda.samples(jags, c('phi'), 20000, thin=10)
auxmcmc <- coda.samples(jags, c('lambda', 'sigma', 'b1','b2','b3', 'theta', 'rho'), 20000, thin=10)
out_jags <- list(phimcmc=phimcmc, auxmcmc=auxmcmc, theta.select=theta.select, rho.select=rho.select, phi.select=phi.select)
save(out_jags,file="./models/m301/out_jags.Rdata")
